#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define M 15

double func(double x,int i)
{
	return pow(x,i);
}
double gunc(double x,double y[],int n)
{
	int i;
	double r = 0;
	for(i = 0;i < n;i++)
		r += y[i]*pow(x,i);
	return r;
}
void fitting(double x[],double b[],int n,int m)
{
	double a[m][n];	//a + bt + ct^2
	double g[n][n];
	double f[n];
	double t[n];
	double y[n];
	int i,j,k;
	//形成A矩阵
	for (i = 0; i < m; i++)
		for (j = 0; j < n; j++)
			a[i][j] = func(x[i],j);
	//形成G矩阵和f向量
	for (i = 0; i < n; i++)
	{
		f[i] = 0;
		for (j = 0; j < n; j++)
		{
			g[i][j] = 0;
			for(k = 0;k < m;k++)
				g[i][j] += a[k][i]*a[k][j];
		}
		for(j = 0;j < m;j++)
			f[i] += a[j][i]*b[j];
	}
	//形成l下三角矩阵
	for(j = 0;j < n;j++)
	{
		for(k = 0;k < j;k++)
			g[j][j] -= g[j][k]*g[j][k];

		g[j][j] = sqrt(g[j][j]);
		for(i = j+1;i < n;i++)
		{
			for(k = 0;k < j;k++)
				g[i][j] -= g[i][k]*g[j][k];
			g[i][j] /= g[j][j];
		}
	}
	//计算L t = f中的t
	for(i = 0;i < n;i++)
	{
		t[i] = f[i];
		for(j = 0;j < i;j++)
			t[i] -= g[i][j]*t[j];
		t[i] /= g[i][i];
	}
	//计算L(转置) y = t中的y
	for(i = n-1;i >= 0;i--)
	{
		y[i] = t[i];
		for(j = n-1;j > i;j--)
			y[i] -= g[j][i]*y[j];
		y[i] /= g[i][i];
	}
	for(i = 0;i < n;i++)
		printf("系数%d：%f\n",i,y[i]);
	//for(i = 0;i < m;i++)
	//	printf("%f %f\n",b[i],gunc(x[i],y,n));
}
int main()
{
	int m = M;
	double x[M] = {1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8};
	double b[M] = {33.4,79.50,122.65,159.05,189.15,214.15,238.65,252.20,267.55,280.50,296.65,301.65,310.40,318.15,325.15};
	//计算y = a + bt + ct^2的系数
	printf("y = a + bt + ct^2\n");
	fitting(x,b,3,m);
	//计算y = ae^bt的系数，转换为ln(y) = ln(a) + bt
	//处理原始数据ln(y)
	double lb[m]; int i;
	for(i = 0;i < m;i++)
		lb[i] = log(b[i]);
	printf("\nln(y) = ln(a) + bt\n");
	fitting(x,lb,2,m);
	return 0;
}
